/* grid.c is part of VIP program
 * Written by J. Salvador Arias, 2010
 *      CONICET, INSUE, Instituto Miguel Lillo,
 *      Miguel Lillo 205, S.M. de Tucuman (4000)
 *      Tucuman, Argentina
 *      e-mail: jsalarias@csnat.unt.edu.ar
 *
 * VIP is a free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public
 * License as published by the Free Software Foundation and
 * Creative Commons.
 *
 * This file and the VIP program is distributed in the hope
 * that it will be useful, but WITHOUT ANY WARRANTY; without
 * even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * To view a copy of this license, visit
 * http://creativecommons.org/licenses/GPL/2.0/
 * or send a letter to
 * Creative Commons, 171 Second Street, Suite 300,
 * San Francisco, California, 94105, USA.
 */

#include "grid.h"

static GNode* qTree = NULL;
static GPtrArray* cellArray = NULL;

static GNode* CreateQTree (GNode* anc, GNode* sibling, SINTEGER iniRow, SINTEGER endRow, SINTEGER iniCol, SINTEGER endCol,
                           INTARRAY* matrix);
static void MakeQTree ();

static void MakeCellObs (SCOUNTER destCell);

static void MapQNodeDown (GNode* qNode, BITFIELD* setA, BITFIELD* setB);
static void MapQNodeUp (GNode* qNode, BITFIELD* setA, BITFIELD* setB);

/* Phylogenetic diversity */
static SCOUNTER PhyloDivNode (GNode* node, PHYLONODE* anc, SCOUNTER cell);

/* Here I use a recursive function to build the qTree, becuase it is more difficult to build a non-recursive option! */
GNode* CreateQTree (GNode* anc, GNode* sibling, SINTEGER iniRow, SINTEGER endRow, SINTEGER iniCol, SINTEGER endCol,
                    INTARRAY* matrix) {
 QUADTREE* qNode;
 GNode* actNode;
 GNode* nuSibling;
 SCOUNTER r, c, actArea;
 SCOUNTER halfRow, halfCol;

 UpdateOnCalculus ();
 qNode = g_slice_alloc (sizeof (QUADTREE));
 qNode -> firstRow = iniRow;
 qNode -> lastRow = endRow;
 qNode -> firstCol = iniCol;
 qNode -> lastCol = endCol;
 qNode -> size = (endRow - iniRow + 1) * (endCol - iniCol + 1);
 qNode -> state = FALSE;
 qNode -> obsCell = FALSE;
 qNode -> cell = -1;
 actNode = g_node_new (qNode);
 AddChildren (actNode, anc, sibling);

 for (r = iniRow; r <= endRow; ++ r) {
   for (c = iniCol; c <= endCol; ++ c) {
     actArea = r * numCols + c;
     if (GetFromArray (matrix, actArea) > 0) {
       qNode -> state = TRUE;
       break;
     }
   }
   if (qNode -> state == TRUE) break;
 }

 if (qNode -> size == 1) {
   if (qNode -> state == TRUE) {
     qNode -> cell = numCells;
     g_ptr_array_add (cellArray, qNode);
     ++ numCells;
   }
   return actNode;
 }
 if (qNode -> state == FALSE) return actNode;
 nuSibling = NULL;
 if ((endCol > iniCol) && (endRow > iniRow)) {
   halfCol = iniCol + ((endCol - iniCol) / 2);
   halfRow = iniRow + ((endRow - iniRow) / 2);
   nuSibling = CreateQTree (actNode, nuSibling, iniRow, halfRow, iniCol, halfCol, matrix);
   nuSibling = CreateQTree (actNode, nuSibling, iniRow, halfRow, halfCol + 1, endCol, matrix);
   nuSibling = CreateQTree (actNode, nuSibling, halfRow + 1, endRow, iniCol, halfCol, matrix);
   CreateQTree (actNode, nuSibling, halfRow + 1, endRow, halfCol + 1, endCol, matrix);
 }
 else if (endCol > iniCol) {
   halfCol = iniCol + ((endCol - iniCol) / 2);
   nuSibling = CreateQTree (actNode, nuSibling, iniRow, endRow, iniCol, halfCol, matrix);
   CreateQTree (actNode, nuSibling, iniRow, endRow, halfCol + 1, endCol, matrix);
 }
 else if (endRow > iniRow) {
   halfRow = iniRow + ((endRow - iniRow) / 2);
   nuSibling = CreateQTree (actNode, nuSibling, iniRow, halfRow, iniCol, endCol, matrix);
   CreateQTree (actNode, nuSibling, halfRow + 1, endRow, iniCol, endCol, matrix);
 }
 return actNode;
}

SCOUNTER GetCell (SINTEGER row, SINTEGER col) {
 GNode* actNode;
 QUADTREE* qNode;

 qNode = NULL;
 for (actNode = qTree; actNode != NULL;) {
   qNode = (QUADTREE*) actNode -> data;
   if ((col >= qNode -> firstCol) && (col <= qNode -> lastCol) && (row >= qNode -> firstRow) && (row <= qNode -> lastRow)) {
     if (qNode -> size == 1) break;
     actNode = g_node_first_child (actNode);
     qNode = NULL;
     continue;
   }
   if (g_node_next_sibling (actNode) == NULL) {
     qNode = NULL;
     break;
   }
   actNode = g_node_next_sibling (actNode);
 }
 if (qNode == NULL) return -1;
 return qNode -> cell;
}

void MakeQTree () {
 SCOUNTER maxAreas, actArea;
 SCOUNTER row, col, r, c, num;
 GEOREF* actRec;
 GList* navGeo;
 INTARRAY* matrix;
 SREAL x, y;

 maxAreas = numCols *  numRows;
 matrix = CreateIntArray (maxAreas);
 num = 0;
 for (actArea = 0; actArea < maxAreas; ++ actArea)
   AddToArray (matrix, num);
 for (navGeo = GETGEOREF; navGeo != NULL; navGeo = g_list_next (navGeo)) {
   actRec = (GEOREF*) navGeo -> data;
   x = actRec -> azimuth - gridStartX;
   y = actRec -> inclination - gridStartY;
   col = (SINTEGER) (x / actualGridSize);
   row = (SINTEGER) (y / actualGridSize);
   for (r = row - actualFill; r <= (row + actualFill); ++ r) {
     if (r < 0) continue;
     if (r >= numRows) break;
     for (c = col - actualFill; c <= (col + actualFill); ++ c) {
       if (c < 0) continue;
       if (c >= numCols) break;
       actArea = (r * numCols) + c;
       num = GetFromArray (matrix, actArea);
       ++ num;
       InsertValToArray (matrix, num, actArea);
     }
   }
 }
 numCells = 0;
 cellArray = g_ptr_array_new ();
 qTree = CreateQTree (NULL, NULL, 0, numRows - 1, 0, numCols - 1, matrix);
 DestroyIntArray (matrix);
}

void CreateGrid () {
 GList* navTerm;
 GList* navGeo;
 TERMINAL* actTerm;
 PHYLONODE* node;
 GEOREF* actGeo;
 SCOUNTER actCell;
 SCOUNTER row, col;
 SREAL x, y;
 SCOUNTER termFill;

 numCols = ((GETMAXAZM - GETMINAZM) / actualGridSize) + 2;
 numRows = ((GETMAXINC - GETMININC) / actualGridSize) + 2;
 gridStartX = ((SINTEGER) (GETMINAZM / actualGridSize)) * actualGridSize;
 gridStartY = ((SINTEGER) (GETMININC / actualGridSize)) * actualGridSize;

 if (qTree == NULL) {
   MakeQTree ();
   StartBitFields (numCells);
   for (navTerm = GETSPP; navTerm != NULL; navTerm = g_list_next (navTerm)) {
     actTerm = (TERMINAL*) navTerm -> data;
     actTerm -> obs = CreateBitField ();
     node = (PHYLONODE*) actTerm -> node -> data;
     termFill = (node -> fill > actualFill) ? actualFill : node -> fill;
     for (navGeo = actTerm -> firstGeo; navGeo != NULL; navGeo = g_list_next (navGeo)) {
       UpdateOnCalculus ();
       actGeo = (GEOREF*) navGeo -> data;
       if (actGeo -> sp != navTerm) break;
       x = actGeo -> azimuth - gridStartX;
       y = actGeo -> inclination - gridStartY;
       col = (SINTEGER) (x / actualGridSize);
       row = (SINTEGER) (y / actualGridSize);
       actCell = GetCell (row, col);
       actGeo -> actCell = actCell;
       if (IsBitOn (actTerm -> obs, actCell) == TRUE) continue;
       MakeCellObs (actCell);
       PutBitOn (actTerm -> obs, actCell);
     }
   }
 }
 MakeOptimization (FALSE);
}

void DestroyGrid () {
 GList* navTerm;
 GList* navGeo;
 GNode* actNode;
 GNode* nuAct;
 TERMINAL* actTerm;
 GEOREF* actGeo;

 if (qTree == NULL) return;
 UpdateOnCalculus ();
 DestroyRecs ();
 RemoveAlloSets ();
 for (navTerm = GETSPP; navTerm != NULL; navTerm = g_list_next (navTerm)) {
   UpdateOnCalculus ();
   actTerm = (TERMINAL*) navTerm -> data;
   if (actTerm -> obs != NULL)
     g_slice_free1 (GETBITFIELDSIZE, actTerm -> obs);
   actTerm -> obs = NULL;
 }
 for (navGeo = GETGEOREF; navGeo != NULL; navGeo = g_list_next (navGeo)) {
   UpdateOnCalculus ();
   actGeo = (GEOREF*) navGeo -> data;
   actGeo -> actCell = -1;
 }
 actNode = qTree;
 while (actNode != NULL) {
   if (g_node_first_child (actNode) != NULL) {
     actNode = g_node_last_child (actNode);
     continue;
   }
   if (g_node_prev_sibling (actNode) != NULL)
     nuAct = g_node_prev_sibling (actNode);
   else
     nuAct = actNode -> parent;
   g_node_unlink (actNode);
   g_slice_free1 (sizeof (QUADTREE), actNode -> data);
   g_node_destroy (actNode);
   actNode = nuAct;
 }
 qTree = NULL;
 g_ptr_array_free (cellArray, TRUE);
 cellArray = NULL;
}

SFLAG GetGridState () {
 return (cellArray != NULL) ? TRUE : FALSE;
}

void GetColRowFromCell (SCOUNTER cell, SINTEGER* row, SINTEGER* col) {
 QUADTREE* qNode;

 qNode = (QUADTREE*) g_ptr_array_index (cellArray, cell);
 *row = qNode -> firstRow;
 *col = qNode -> firstCol;
}

void MapIntoQTree (BITFIELD* setA, BITFIELD* setB) {
 MapQNodeDown (qTree, setA, setB);
 MapQNodeUp (qTree, setA, setB);
}

void MapQNodeDown (GNode* qNode, BITFIELD* setA, BITFIELD* setB) {
 QUADTREE* quad;
 QUADTREE* child;
 GNode* navChild;

 quad = (QUADTREE*) qNode -> data;
 quad -> state = 0;
 if (quad -> size == 1) {
   if (quad -> cell != -1) {
     if (IsBitOn (setA, quad -> cell))
       quad -> state |= 1;
     if (IsBitOn (setB, quad -> cell))
       quad -> state |= 2;
   }
   return;
 }
 for (navChild = g_node_first_child (qNode); navChild != NULL; navChild = g_node_next_sibling (navChild)) {
   MapQNodeDown (navChild, setA, setB);
   child = (QUADTREE*) navChild -> data;
   quad -> state |= child -> state;
 }
 return;
}

void MapQNodeUp (GNode* qNode, BITFIELD* setA, BITFIELD* setB) {
 QUADTREE* quad;
 QUADTREE* anc;
 GNode* navChild;

 quad = (QUADTREE*) qNode -> data;

 if ((quad -> state == 0) && (qNode -> parent != NULL)) {
   anc = (QUADTREE*) qNode -> parent -> data;
   quad -> state = anc -> state;
 }
 for (navChild = g_node_first_child (qNode); navChild != NULL; navChild = g_node_next_sibling (navChild))
   MapQNodeUp (navChild, setA, setB);
}

SFLAG GetRowColState (SINTEGER row, SINTEGER col) {
 GNode* actNode;
 QUADTREE* qNode;

 qNode = NULL;
 for (actNode = qTree; actNode != NULL;) {
   qNode = (QUADTREE*) actNode -> data;
   if ((col >= qNode -> firstCol) && (col <= qNode -> lastCol) && (row >= qNode -> firstRow) && (row <= qNode -> lastRow)) {
     if (qNode -> state != 3) break;
     if (g_node_first_child (actNode) == NULL) break;
     actNode = g_node_first_child (actNode);
     qNode = NULL;
     continue;
   }
   if (g_node_next_sibling (actNode) == NULL) {
     qNode = NULL;
     break;
   }
   actNode = g_node_next_sibling (actNode);
 }
 if (qNode == NULL) return 0;
 return qNode -> state;
}

void MakeCellObs (SCOUNTER destCell) {
 QUADTREE* qNode;

 qNode = (QUADTREE*) g_ptr_array_index (cellArray, destCell);
 qNode -> obsCell = TRUE;
}

SCOUNTER GetCellIfObs (SINTEGER row, SINTEGER col) {
 GNode* actNode;
 QUADTREE* qNode;

 qNode = NULL;
 for (actNode = qTree; actNode != NULL;) {
   qNode = (QUADTREE*) actNode -> data;
   if ((col >= qNode -> firstCol) && (col <= qNode -> lastCol) && (row >= qNode -> firstRow) && (row <= qNode -> lastRow)) {
     if (qNode -> size == 1) break;
     actNode = g_node_first_child (actNode);
     qNode = NULL;
     continue;
   }
   if (g_node_next_sibling (actNode) == NULL) {
     qNode = NULL;
     break;
   }
   actNode = g_node_next_sibling (actNode);
 }
 if (qNode == NULL) return -1;
 return (qNode -> obsCell == TRUE) ? qNode -> cell : -1;
}

void GetObsSet (BITFIELD* obs, GNode* refNode) {
 GNode* actNode;
 PHYLONODE* node;
 PHYLONODE* anc;
 TERMINAL* term;

 node = (PHYLONODE*) refNode -> data;
 ResetBitField (obs);
 for (actNode = refNode; actNode != NULL;) {
   if (IsDesc (actNode, refNode) != TRUE) break;
   node = (PHYLONODE*) actNode -> data;
   if (actNode -> parent != NULL)
     anc = (PHYLONODE*) actNode -> parent -> data;
   if ((actNode != refNode) && ((node -> rec -> numCells == 0) ||
                 (node -> rec -> nodeFlag & REMOV_NODE) || (node == anc -> rec -> removSet))) {
     /* A removed distribution, ignores all of its descendants */
     while (actNode != NULL) {
       if (g_node_next_sibling (actNode) != NULL) {
         actNode = g_node_next_sibling (actNode);
         break;
       }
       actNode = actNode -> parent;
     }
     continue;
   }
   if (g_node_first_child (actNode) != NULL) {
   /* Not a term, goes to the first descendant */
     actNode = g_node_first_child (actNode);
     continue;
   }
   term = node -> term;
   BitFieldORsum (obs, term -> obs);
   while (actNode != NULL) {
     if (g_node_next_sibling (actNode) != NULL) {
       actNode = g_node_next_sibling (actNode);
       break;
     }
     actNode = actNode -> parent;
   }
 }
}

void DoFill (BITFIELD* dest, BITFIELD* obs, SCOUNTER fillValue) {
 SCOUNTER navCell;
 SCOUNTER actCell;
 SCOUNTER row, col, fill;
 SINTEGER r, c;
 SINTEGER xBlocks, yBlocks, blocks;

 ResetBitField (dest);
 fill = (fillValue < actualFill) ? fillValue : actualFill;
 for (navCell = 0; navCell <= numCells; ++ navCell) {
   if (IsBitOn (obs, navCell) != TRUE) continue;
   GetColRowFromCell (navCell, &row, &col);
   for (r = row - fill; r <= (row + fill); ++ r) {
     if (r < 0) continue;
     if (r >= numRows) break;
     for (c = col - fill; c <= (col + fill); ++ c) {
       if (c < 0) continue;
       if (c >= numCols) break;
       if (gridMode & VONNEUMANN_NEIGHBORS) {
         xBlocks = (r < row) ? row - r : r - row;
         yBlocks = (c < col) ? col - c : c - col;
         blocks = xBlocks + yBlocks;
         if (blocks > fill) continue;
       }
       actCell = GetCell (r, c);
       PutBitOn (dest, actCell);
     }
   }
 }
}

SCOUNTER PhyloDivNode (GNode* node, PHYLONODE* anc, SCOUNTER cell) {
 PHYLONODE* actNode;
 GNode* navDesc;
 SCOUNTER toRet;

 UpdateOnCalculus ();
 actNode = (PHYLONODE*) node -> data;
 if (IsBitOn (actNode -> rec -> cells, cell) != TRUE) return 0;
 if (anc != NULL)
   toRet = anc -> level - actNode -> level;
 else
   toRet = 0;
 if (g_node_first_child (node) != NULL) {
   for (navDesc = g_node_first_child (node); navDesc != NULL; navDesc = g_node_next_sibling (navDesc))
     toRet += PhyloDivNode (navDesc, actNode, cell);
 }
 return toRet;
}

SCOUNTER CalculatePhyloDiv (SCOUNTER cell, SCOUNTER max) {
 SCOUNTER toRet, thisTree;
 SCOUNTER extra;
 GList* navTree;
 PHYLOTREE* actTree;
 PHYLONODE* root;

 if (cellArray == NULL) return 0;
 toRet = 0;
 for (navTree = GETTREES; navTree != NULL; navTree = g_list_next (navTree)) {
   actTree = (PHYLOTREE*) navTree -> data;
   thisTree = PhyloDivNode (actTree -> root, NULL, cell);
   if (thisTree > 0) {
     root = (PHYLONODE*) actTree -> root -> data;
     extra = max - root -> level;
     thisTree += extra;
   }
   toRet += thisTree;
 }
 return toRet;
}
